library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(knitr)
library(ggplot2)
library(dplyr)
library(lubridate)
library(MMWRweek)
library(gridExtra)
##
## Attaching package: 'gridExtra'
##
## The following object is masked from 'package:dplyr':
##
## combine
library(readxl)
library(excessmort)
library(plotly)
##
## Attaching package: 'plotly'
##
## The following object is masked from 'package:ggplot2':
##
## last_plot
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following object is masked from 'package:graphics':
##
## layout
library(patchwork)
library(scales)
##
## Attaching package: 'scales'
##
## The following object is masked from 'package:purrr':
##
## discard
##
## The following object is masked from 'package:readr':
##
## col_factor
library(cowplot)
##
## Attaching package: 'cowplot'
##
## The following object is masked from 'package:patchwork':
##
## align_plots
##
## The following object is masked from 'package:lubridate':
##
## stamp
library(broom)
library(gt)
## Warning: Can't find generic `as.gtable` in package gtable to register S3 method.
## ℹ This message is only shown to developers using devtools.
## ℹ Do you need to update gtable to the latest version?
##
## Attaching package: 'gt'
##
## The following object is masked from 'package:cowplot':
##
## as_gtable
library(webshot2)
library(stringr)
library(purrr)
library(kableExtra)
##
## Attaching package: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rows
source("~/Desktop/Thesis Work/Hospitalization/11-26 data/Deficit_Model_1.0.R")
allcause_state_byweek <- read_excel("allcause_state_byweek_suppressed_25mar21.xlsx")
icd10categories <- read_excel("icd10categories_updated_321.xlsx")
#weekly state data
all_cause_state_byweek <- allcause_state_byweek %>%
rename(weekly_count_statelevel = count)
all_cause_state_by_week <- all_cause_state_byweek %>%
separate(mmwr_week, into = c("MMWRyear", "MMWRweek"), sep = "-", convert = TRUE) %>%
mutate(
date = MMWRweek2Date(MMWRyear = MMWRyear, MMWRweek = MMWRweek, MMWRday = 1),
statecount_numeric = as.numeric(ifelse(weekly_count_statelevel == "<5", NA, weekly_count_statelevel)),
statecount_suppressed = weekly_count_statelevel == "<5",
category = as.character(category)
)
all_cause_state_by_week %>%
group_by(category) %>%
reframe(missing_total = sum(statecount_suppressed),
count_total = sum(statecount_numeric, na.rm = TRUE))
## # A tibble: 10 × 3
## category missing_total count_total
## <chr> <int> <dbl>
## 1 cardiovascular 0 373147
## 2 covid 1 58262
## 3 digestive 0 182701
## 4 infectious 0 112279
## 5 infectious respiratory 0 29816
## 6 injury 0 331205
## 7 other 0 1740863
## 8 pneumonia 0 97954
## 9 renal 0 102199
## 10 respiratory 0 164424
all_cause_state_by_week
## # A tibble: 3,140 × 7
## MMWRyear MMWRweek category weekly_count_statelevel date
## <int> <int> <chr> <chr> <date>
## 1 2019 1 cardiovascular 843 2018-12-30
## 2 2019 1 digestive 441 2018-12-30
## 3 2019 1 infectious 299 2018-12-30
## 4 2019 1 infectious respiratory 240 2018-12-30
## 5 2019 1 injury 691 2018-12-30
## 6 2019 1 other 4090 2018-12-30
## 7 2019 1 pneumonia 427 2018-12-30
## 8 2019 1 renal 220 2018-12-30
## 9 2019 1 respiratory 522 2018-12-30
## 10 2019 1 covid 0 2018-12-30
## # ℹ 3,130 more rows
## # ℹ 2 more variables: statecount_numeric <dbl>, statecount_suppressed <lgl>
str(icd10categories)
## tibble [869 × 5] (S3: tbl_df/tbl/data.frame)
## $ icd10code : chr [1:869] "A00" "A01" "A02" "A03" ...
## $ category_original : chr [1:869] "infectious" "infectious" "infectious" "infectious" ...
## $ category_1113update: chr [1:869] "infectious" "infectious" "infectious" "infectious" ...
## $ category_123update : chr [1:869] "infectious" "infectious" "infectious" "infectious" ...
## $ category_0321update: chr [1:869] "infectious" "infectious" "infectious" "infectious" ...
names(icd10categories)
## [1] "icd10code" "category_original" "category_1113update"
## [4] "category_123update" "category_0321update"
dim(icd10categories)
## [1] 869 5
summary(icd10categories)
## icd10code category_original category_1113update category_123update
## Length:869 Length:869 Length:869 Length:869
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
## category_0321update
## Length:869
## Class :character
## Mode :character
sapply(icd10categories, class)
## icd10code category_original category_1113update category_123update
## "character" "character" "character" "character"
## category_0321update
## "character"
glimpse(icd10categories)
## Rows: 869
## Columns: 5
## $ icd10code <chr> "A00", "A01", "A02", "A03", "A04", "A05", "A06", "…
## $ category_original <chr> "infectious", "infectious", "infectious", "infecti…
## $ category_1113update <chr> "infectious", "infectious", "infectious", "infecti…
## $ category_123update <chr> "infectious", "infectious", "infectious", "infecti…
## $ category_0321update <chr> "infectious", "infectious", "infectious", "infecti…
colSums(is.na(icd10categories))
## icd10code category_original category_1113update category_123update
## 0 418 418 1
## category_0321update
## 0
str(allcause_state_byweek)
## tibble [3,140 × 3] (S3: tbl_df/tbl/data.frame)
## $ mmwr_week: chr [1:3140] "2019-01" "2019-01" "2019-01" "2019-01" ...
## $ category : chr [1:3140] "cardiovascular" "digestive" "infectious" "infectious respiratory" ...
## $ count : chr [1:3140] "843" "441" "299" "240" ...
names(allcause_state_byweek)
## [1] "mmwr_week" "category" "count"
dim(allcause_state_byweek)
## [1] 3140 3
summary(allcause_state_byweek)
## mmwr_week category count
## Length:3140 Length:3140 Length:3140
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
sapply(allcause_state_byweek,class)
## mmwr_week category count
## "character" "character" "character"
glimpse(allcause_state_byweek)
## Rows: 3,140
## Columns: 3
## $ mmwr_week <chr> "2019-01", "2019-01", "2019-01", "2019-01", "2019-01", "2019…
## $ category <chr> "cardiovascular", "digestive", "infectious", "infectious res…
## $ count <chr> "843", "441", "299", "240", "691", "4090", "427", "220", "52…
colSums(is.na(allcause_state_byweek))
## mmwr_week category count
## 0 0 0
deficit_plot <- function(fit,
title = "",
ylim = NULL,
show.data = TRUE,
alpha = 0.05){
requireNamespace("ggplot2")
z <- qnorm(1 - alpha/2)
p <- with(fit, data.frame(date = date,
y = (observed - expected)/expected,
increase = fitted,
sd = sd,
se = se)) %>%
ggplot(aes(date, y)) +
geom_ribbon(aes(ymin = increase - z * se, ymax = increase + z * se), alpha = 0.5, fill = "#3366FF") + #needs to be edited
geom_hline(yintercept = 0) +
ylab("% Change from expected count") +
scale_y_continuous(labels = scales::percent) +
scale_x_date(date_labels = "%b %Y") +
ggtitle(title) +
xlab("Date")
if(show.data) p <- p + geom_point(alpha = 0.3)
if(!is.null(ylim)) p <- p + coord_cartesian(ylim = ylim)
return(p + geom_line(aes(y = increase), size=0.70, col="#3366FF"))
}
deficit_cumulative <- function(fit, start, end){
if (!"curve_fit" %in% attr(fit, "type"))
stop("This is not the correct deficit_model fit, needs curve fit.")
ind <- which(fit$date %in% seq(start, end, by = "day"))
n <- length(ind) # returns the # of elements in the ind object
A <- matrix(1, n, n)
A[upper.tri(A)] <- 0 # It sets all the upper triangular elements of matrix A (excluding the diagonal) to 0
A <- sweep(A, 2, fit$expected[ind], FUN = "*") # This line multiplies each column (2 means columns) of the matrix A by a corresponding value from fit$expected[ind]
# So, now it's scaling each column of A by expected values from the model fit (fit$expected[ind])
fhat <- matrix(fit$fitted[ind], ncol = 1) # creates a column vector (matrix with 1 column) called fhat from the values in fit$fitted[ind]
# Below, this code is constructing a data frame with:
fit_deficit <- A %*% fhat # multiply design matrix A by fitted values fhat to get cumulative modeled effect.
obs_deficit <- cumsum(fit$observed[ind] - fit$expected[ind])
varcov <- fit$x[ind,] %*% fit$betacov %*% t(fit$x[ind,]) # Compute variance-covariance matrix of the linear predictor based on model design matrix x and estimated betacov
diag(varcov) <- fit$se[ind]^2 # Replace diagonal entries of varcov with squared standard errors — probably to correct for over/under-dispersion or model misspecification.
fitted_se <- sqrt(diag(A %*% varcov %*% t(A)))
sd <- sqrt(diag(A %*% (fit$cov[ind,ind] + diag(fit$log_expected_se[ind]^2)) %*% t(A))) # this computes the variance of the observed cumulative deficit.
# It adds fit$cov (covariance of observed data) and extra uncertainty from fit$log_expected_se.
data.frame(date = fit$date[ind], # The time points
observed = obs_deficit, # Cumulative observed - expected difference
sd = sd, # Modeled cumulative deficit/excess from the model.
fitted = fit_deficit, # Standard deviation of the observed cumulative deficit.
se = fitted_se) # Standard error of the fitted cumulative deficit.
}
exclude_dates <- seq(lubridate::make_date(2020, 3, 1),
lubridate::make_date(2021, 2, 20),
by = "day")
# cardiovascular data
cardio <- all_cause_state_by_week %>%
filter(category == "cardiovascular") %>%
filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
outcome = statecount_numeric)
3
## [1] 3
# deficit model with COVID-19 dates excluded
cardio_deficit <- compute_expected(cardio,
exclude = exclude_dates,
frequency = 52)
## Overall death rate is 8.85.
cardio_deficit
## # A tibble: 312 × 12
## MMWRyear MMWRweek category weekly_count_statelevel date
## <int> <int> <chr> <chr> <date>
## 1 2019 2 cardiovascular 1233 2019-01-06
## 2 2019 3 cardiovascular 1160 2019-01-13
## 3 2019 4 cardiovascular 1201 2019-01-20
## 4 2019 5 cardiovascular 1157 2019-01-27
## 5 2019 6 cardiovascular 1286 2019-02-03
## 6 2019 7 cardiovascular 1261 2019-02-10
## 7 2019 8 cardiovascular 1141 2019-02-17
## 8 2019 9 cardiovascular 1156 2019-02-24
## 9 2019 10 cardiovascular 1208 2019-03-03
## 10 2019 11 cardiovascular 1232 2019-03-10
## # ℹ 302 more rows
## # ℹ 7 more variables: statecount_numeric <dbl>, statecount_suppressed <lgl>,
## # population <dbl>, outcome <dbl>, log_expected_se <dbl>, expected <dbl>,
## # excluded <lgl>
cardio_deficit <- cardio_deficit %>%
mutate(date = as.Date(paste(MMWRyear, MMWRweek, 1, sep = "-"), "%Y-%U-%u"))
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `date = as.Date(paste(MMWRyear, MMWRweek, 1, sep = "-"),
## "%Y-%U-%u")`.
## Caused by warning in `strptime()`:
## ! (0-based) yday 369 in year 2020 is invalid
combined_plot_1 <- ggplot(cardio_deficit, aes(x = date)) +
geom_line(aes(y = expected), color = "blue", size = 1.5) + # Expected outcome line
geom_point(aes(y = outcome, color = factor(excluded)), size = 2) + # Observed outcomes
scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
labels = c("Included in Model", "Excluded from Model")) +
labs(title = "Weekly Cardiovascular Hospitalizations in MA",
x = "Date",
y = "Weekly Hospitilizations",
color = "Data Point Status") +
theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
combined_plot_1
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
ggsave("combined_plot_1.png", plot = combined_plot_1, width = 8, height = 6, dpi = 300)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
fit_x <- cardio_deficit %>%
deficit_model(exclude = exclude_dates,
start = as.Date("2019-01-14"),
end = as.Date("2024-12-31"),
knots.per.year = 12,
verbose = FALSE)
# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_x <- deficit_plot(fit_x,
title = "Deviations from Expected Counts for Cardio in MA"
)
deficit_plot_x
fit_x$detected_intervals
## start end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-09-23 2019-12-02 8.144139 8.669496 0.07650853 12062
## 2 2020-01-20 2020-02-03 8.276251 8.703193 0.14678711 3343
## 3 2020-03-02 2021-02-15 7.265876 8.874462 0.03594969 49893
## 4 2021-12-27 2022-01-31 8.304721 9.024090 0.10569035 6709
## 5 2023-05-01 2023-05-01 9.677494 9.573722 0.26665497 1303
## expected deficit sd
## 1 12840.088 -778.0880 113.31411
## 2 3515.453 -172.4533 59.29126
## 3 60938.774 -11045.7736 246.85780
## 4 7290.145 -581.1446 85.38234
## 5 1289.028 13.9721 35.90303
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit <- deficit_cumulative(fit_x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_1 <- cumulative_deficit %>%
ggplot(aes(date)) +
geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
geom_line(aes(y = observed),
color = "white",
size = 1) +
geom_line(aes(y = observed)) +
geom_point(aes(y = observed)) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Date",
y = "Cumulative Deficit Hospitilization-Cardio",
title = "Cumulative Deficit Hospitilization for Cardio Cases in MA",
subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_1
library(patchwork)
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")
combined_plot_1_adj <- combined_plot_1 +
theme(legend.position = "bottom",
legend.justification = c(1, 0),
legend.box.just = "right") +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
deficit_plot_x_adj <- deficit_plot_x +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_1_adj <- cumulative_deficit_plot_1 +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
final_combined_plot_1 <- combined_plot_1_adj /
deficit_plot_x_adj /
cumulative_deficit_plot_1_adj +
plot_layout(ncol = 1, heights = c(1, 1, 1))
final_combined_plot_1
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
ggsave("Cardio_combined_figure.png",
plot = final_combined_plot_1,
width = 11,
height = 10.5,
dpi = 300,
bg = "white")
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
ggsave("cardio_combined_figure.png",
plot = final_combined_plot_1,
width = 11, # good for standard document width
height = 10.5, # leaves space for three stacked plots
dpi = 300,
bg = "white")
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1),
lubridate::make_date(2021, 2, 27),
by = "day")
# digestive data
digestive <- all_cause_state_by_week %>%
filter(category == "digestive") %>%
filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
outcome = statecount_numeric)
# excess model with COVID-19 dates excluded
digestive_deficit <- compute_expected(digestive,
exclude = exclude_dates,
frequency = 52)
## Overall death rate is 4.33.
digestive_deficit
## # A tibble: 312 × 12
## MMWRyear MMWRweek category weekly_count_statelevel date
## <int> <int> <chr> <chr> <date>
## 1 2019 2 digestive 535 2019-01-06
## 2 2019 3 digestive 570 2019-01-13
## 3 2019 4 digestive 553 2019-01-20
## 4 2019 5 digestive 576 2019-01-27
## 5 2019 6 digestive 609 2019-02-03
## 6 2019 7 digestive 599 2019-02-10
## 7 2019 8 digestive 602 2019-02-17
## 8 2019 9 digestive 562 2019-02-24
## 9 2019 10 digestive 550 2019-03-03
## 10 2019 11 digestive 622 2019-03-10
## # ℹ 302 more rows
## # ℹ 7 more variables: statecount_numeric <dbl>, statecount_suppressed <lgl>,
## # population <dbl>, outcome <dbl>, log_expected_se <dbl>, expected <dbl>,
## # excluded <lgl>
# Plotting
combined_plot_2 <- ggplot(digestive_deficit, aes(x = date)) +
geom_line(aes(y = expected), color = "blue", size = 1.5) + # Expected outcome line
geom_point(aes(y = outcome, color = factor(excluded)), size = 2) + # Observed outcomes
scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
labels = c("Included in Model", "Excluded from Model")) +
labs(title = "Weekly digestive ED Hospitilizations in MA",
x = "Date",
y = "Weekly Hospitilizations",
color = "Data Point Status") +
theme_minimal()
combined_plot_2
fit_2x <- digestive_deficit %>%
deficit_model(exclude = exclude_dates,
start = min(.$date),
end = max(.$date),
knots.per.year = 12,
verbose = FALSE)
# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_2x <- deficit_plot(fit_2x,
title = "Deviations from Expected Counts for Digestive in MA"
)
deficit_plot_2x
fit_2x$detected_intervals
## start end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-09-15 2019-10-13 3.986860 4.261767 0.07956446 2684
## 2 2020-03-01 2020-09-06 3.519113 4.442692 0.03432838 13267
## 3 2020-10-11 2021-03-14 3.714835 4.193817 0.03680019 11504
## 4 2021-12-05 2022-01-30 3.759757 4.186787 0.05877984 4556
## 5 2022-05-01 2022-05-22 4.368984 4.646495 0.09288423 2353
## 6 2022-09-18 2022-10-02 4.263148 4.543836 0.10606203 1722
## expected deficit sd
## 1 2869.071 -185.0706 53.56371
## 2 16748.877 -3481.8769 129.41745
## 3 12987.299 -1483.2986 113.96183
## 4 5073.467 -517.4671 71.22827
## 5 2502.459 -149.4588 50.02458
## 6 1835.377 -113.3774 42.84131
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit <- deficit_cumulative(fit_2x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_2 <- cumulative_deficit %>%
ggplot(aes(date)) +
geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
geom_line(aes(y = observed),
color = "white",
size = 1) +
geom_line(aes(y = observed)) +
geom_point(aes(y = observed)) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Date",
y = "Cumulative Deficit Hospitilization-Digestive",
title = "Cumulative Deficit Hospitilization for Digestive Cases in MA",
subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_2
library(patchwork)
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")
combined_plot_2_adj <- combined_plot_2 +
theme(legend.position = "bottom",
legend.justification = c(1, 0),
legend.box.just = "right") +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
deficit_plot_2x_adj <- deficit_plot_2x +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_2_adj <- cumulative_deficit_plot_2 +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
final_combined_plot_2 <- combined_plot_2_adj /
deficit_plot_2x_adj /
cumulative_deficit_plot_2_adj +
plot_layout(ncol = 1, heights = c(1, 1, 1))
final_combined_plot_2
ggsave("Digestive_combined_figure.png",
plot = final_combined_plot_2,
width = 11,
height = 10.5,
dpi = 300,
bg = "white")
# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1),
lubridate::make_date(2021, 6, 5),
by = "day")
# infectious data
infectious <- all_cause_state_by_week %>%
filter(category == "infectious") %>%
filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
outcome = statecount_numeric)
# deficit model with COVID-19 dates excluded
infectious_deficit <- compute_expected(infectious,
exclude = exclude_dates,
frequency = 52)
## Overall death rate is 2.66.
# Plotting
combined_plot_3 <- ggplot(infectious_deficit, aes(x = date)) +
geom_line(aes(y = expected), color = "blue", size = 1.5) + # Expected outcome line
geom_point(aes(y = outcome, color = factor(excluded)), size = 2) + # Observed outcomes
scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
labels = c("Included in Model", "Excluded from Model")) +
labs(title = "Weekly Infectious ED Hospitilizations in MA",
x = "Date",
y = "Weekly Hospitilizations",
color = "Data Point Status") +
theme_minimal()
combined_plot_3
fit_3x <- infectious_deficit %>%
deficit_model(exclude = exclude_dates,
start = min(.$date),
end = max(.$date),
knots.per.year = 12,
verbose = FALSE)
deficit_plot_3x <- deficit_plot(fit_3x,
title = "Deviations from Expected Counts for Infectious in MA"
)
deficit_plot_3x
fit_3x$detected_intervals
## start end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2020-05-03 2020-08-02 2.384625 2.718529 0.03797626 4495
## 2 2020-10-04 2021-05-30 2.231734 2.726452 0.02405327 10517
## 3 2022-01-23 2022-02-27 2.530161 2.785967 0.05872481 2044
## 4 2022-12-25 2023-01-22 2.565316 2.820316 0.06472516 1727
## 5 2023-06-11 2023-07-02 2.519639 2.766860 0.07167585 1357
## 6 2024-10-20 2024-12-15 2.525210 2.768028 0.04779399 3060
## expected deficit sd
## 1 5124.405 -629.4047 71.58495
## 2 12848.349 -2331.3491 113.35056
## 3 2250.654 -206.6542 47.44106
## 4 1898.669 -171.6689 43.57372
## 5 1490.146 -133.1456 38.60240
## 6 3354.243 -294.2430 57.91583
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit <- deficit_cumulative(fit_3x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_3 <- cumulative_deficit %>%
ggplot(aes(date)) +
geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
geom_line(aes(y = observed),
color = "white",
size = 1) +
geom_line(aes(y = observed)) +
geom_point(aes(y = observed)) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Date",
y = "Cumulative Deficit Hospitilization-Infectious",
title = "Cumulative Deficit Hospitilization for Infectious Cases in MA",
subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_3
library(patchwork)
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")
combined_plot_3_adj <- combined_plot_3 +
theme(legend.position = "bottom",
legend.justification = c(1, 0),
legend.box.just = "right") +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
deficit_plot_3x_adj <- deficit_plot_3x +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_3_adj <- cumulative_deficit_plot_3 +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
final_combined_plot_3 <- combined_plot_3_adj /
deficit_plot_3x_adj /
cumulative_deficit_plot_3_adj +
plot_layout(ncol = 1, heights = c(1, 1, 1))
final_combined_plot_3
ggsave("Infectious_combined_figure.png",
plot = final_combined_plot_3,
width = 11,
height = 10.5,
dpi = 300,
bg = "white")
# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1),
lubridate::make_date(2021, 6, 12),
by = "day")
# infectious respiratory data
infectious_respiratory <- all_cause_state_by_week %>%
filter(category == "infectious respiratory") %>%
filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
outcome = statecount_numeric)
# excess model with COVID-19 dates excluded
infectious_respiratory_deficit <- compute_expected(infectious_respiratory,
exclude = exclude_dates,
frequency = 52)
## Overall death rate is 0.699.
# Plotting
combined_plot_4 <- ggplot(infectious_respiratory_deficit, aes(x = date)) +
geom_line(aes(y = expected), color = "blue", size = 1.5) + # Expected outcome line
geom_point(aes(y = outcome, color = factor(excluded)), size = 2) + # Observed outcomes
scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
labels = c("Included in Model", "Excluded from Model")) +
labs(title = "Weekly Infectious Respiratory ED Hospitilizations in MA",
x = "Date",
y = "Weekly Hospitilizations",
color = "Data Point Status") +
theme_minimal()
combined_plot_4
# -- Fitting excess model
fit_4x <- infectious_respiratory_deficit %>%
deficit_model(exclude = exclude_dates,
start = min(.$date),
end = max(.$date),
knots.per.year = 12,
verbose = FALSE)
# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_4x <- deficit_plot(fit_4x,
title = "Deviations from Expected Counts for Infectious Respiratory in MA"
)
deficit_plot_4x
# -- Intervals of inordinate mortality found by the excess model
fit_4x$detected_intervals
## start end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-10-13 2019-12-01 0.6340876 1.2580696 0.03417564 683
## 2 2020-03-29 2021-05-30 0.1241042 0.8568669 0.01013141 1036
## 3 2021-10-31 2022-03-20 0.4714432 1.3271874 0.02166535 1333
## 4 2023-01-15 2023-02-26 0.6535837 1.2568241 0.03651720 616
## 5 2023-04-09 2023-04-23 0.4208683 0.7014588 0.04167251 170
## 6 2024-11-03 2024-11-10 0.6424430 1.0037101 0.06105179 173
## expected deficit sd
## 1 1355.1149 -672.11494 36.81189
## 2 7152.9723 -6116.97228 84.57525
## 3 3752.6065 -2419.60646 61.25852
## 4 1184.5517 -568.55175 34.41732
## 5 283.3380 -113.33805 16.83265
## 6 270.2837 -97.28366 16.44031
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit <- deficit_cumulative(fit_4x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
cumulative_deficit_plot_4 <- cumulative_deficit %>%
ggplot(aes(date)) +
geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
geom_line(aes(y = observed),
color = "white",
size = 1) +
geom_line(aes(y = observed)) +
geom_point(aes(y = observed)) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Date",
y = "Cumulative Deficit Hospitilization-infectious respiratory",
title = "Cumulative Deficit Hospitilization for infectious respiratory Cases in MA",
subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_4
library(patchwork)
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")
combined_plot_4_adj <- combined_plot_4 +
theme(legend.position = "bottom",
legend.justification = c(1, 0),
legend.box.just = "right") +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
deficit_plot_4x_adj <- deficit_plot_4x +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_4_adj <- cumulative_deficit_plot_4 +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
final_combined_plot_4 <- combined_plot_4_adj /
deficit_plot_4x_adj /
cumulative_deficit_plot_4_adj +
plot_layout(ncol = 1, heights = c(1, 1, 1))
final_combined_plot_4
ggsave("Infectious Respiratory_combined_figure.png",
plot = final_combined_plot_4,
width = 11, # good for standard document width
height = 10.5, # leaves space for three stacked plots
dpi = 300,
bg = "white")
# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1),
lubridate::make_date(2021, 4, 10),
by = "day")
# injury data
injury <- all_cause_state_by_week %>%
filter(category == "injury") %>%
filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
outcome = statecount_numeric)
# excess model with COVID-19 dates excluded
injury_deficit <- compute_expected(injury,
exclude = exclude_dates,
frequency = 52)
## Overall death rate is 7.86.
# Plotting
combined_plot_5 <- ggplot(injury_deficit, aes(x = date)) +
geom_line(aes(y = expected), color = "blue", size = 1.5) + # Expected outcome line
geom_point(aes(y = outcome, color = factor(excluded)), size = 2) + # Observed outcomes
scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
labels = c("Included in Model", "Excluded from Model")) +
labs(title = "Weekly Injury ED Hospitilizations in MA",
x = "Date",
y = "Weekly Hospitilizations",
color = "Data Point Status") +
theme_minimal()
combined_plot_5
# -- Fitting excess model
fit_5x <- injury_deficit %>%
deficit_model(exclude = exclude_dates,
start = min(.$date),
end = max(.$date),
knots.per.year = 12,
verbose = FALSE)
# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_5x <- excess_plot(fit_5x,
title = "Deviations from Expected Counts for Injuries in MA"
)
deficit_plot_5x
# -- Intervals of inordinate mortality found by the excess model
fit_5x$detected_intervals
## start end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-10-20 2019-11-24 6.756174 7.279676 0.09492698 5458
## 2 2020-01-12 2020-09-13 6.349953 7.623419 0.03965819 30779
## 3 2020-10-04 2021-02-07 6.578836 7.492868 0.05411984 16830
## 4 2021-03-07 2021-04-11 7.120101 7.580370 0.09686766 5752
## 5 2021-05-02 2021-05-16 7.600386 7.944070 0.14023944 3070
## 6 2021-07-25 2021-08-08 7.979167 8.372155 0.14396843 3223
## 7 2021-12-19 2022-01-16 7.140402 7.682093 0.10682282 4807
## 8 2022-07-31 2022-08-21 8.268205 8.668967 0.12687117 4453
## expected deficit sd
## 1 5880.914 -422.9136 76.68712
## 2 36951.644 -6172.6438 192.22810
## 3 19168.280 -2338.2799 138.44956
## 4 6123.830 -371.8301 78.25490
## 5 3208.823 -138.8232 56.64648
## 6 3381.738 -158.7384 58.15272
## 7 5171.673 -364.6727 71.91434
## 8 4668.838 -215.8381 68.32890
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit <- deficit_cumulative(fit_5x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_5 <- cumulative_deficit %>%
ggplot(aes(date)) +
geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
geom_line(aes(y = observed),
color = "white",
size = 1) +
geom_line(aes(y = observed)) +
geom_point(aes(y = observed)) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Date",
y = "Cumulative Deficit Hospitilization-Injury",
title = "Cumulative Deficit Hospitilization for Injury Cases in MA",
subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_5
library(patchwork)
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")
combined_plot_5_adj <- combined_plot_5 +
theme(legend.position = "bottom",
legend.justification = c(1, 0),
legend.box.just = "right") +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
deficit_plot_5x_adj <- deficit_plot_5x +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_5_adj <- cumulative_deficit_plot_5 +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
final_combined_plot_5 <- combined_plot_5_adj /
deficit_plot_5x_adj /
cumulative_deficit_plot_5_adj +
plot_layout(ncol = 1, heights = c(1, 1, 1))
final_combined_plot_5
ggsave("Injury_combined_figure.png",
plot = final_combined_plot_5,
width = 11, # good for standard document width
height = 10.5, # leaves space for three stacked plots
dpi = 300,
bg = "white")
# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1),
lubridate::make_date(2022, 5, 7),
by = "day")
# pneumonia data
pneumonia <- all_cause_state_by_week %>%
filter(category == "pneumonia") %>%
filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
outcome = statecount_numeric)
# excess model with COVID-19 dates excluded
pneumonia_deficit <- compute_expected(pneumonia,
exclude = exclude_dates,
frequency = 52)
## Overall death rate is 2.32.
# Plotting
combined_plot_6 <- ggplot(pneumonia_deficit, aes(x = date)) +
geom_line(aes(y = expected), color = "blue", size = 1.5) + # Expected outcome line
geom_point(aes(y = outcome, color = factor(excluded)), size = 2) + # Observed outcomes
scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
labels = c("Included in Model", "Excluded from Model")) +
labs(title = "Weekly Pneumonia ED Hospitilizations in MA",
x = "Date",
y = "Weekly Hospitilizations",
color = "Data Point Status") +
theme_minimal()
combined_plot_6
# -- Fitting excess model
fit_6x <- pneumonia_deficit %>%
deficit_model(exclude = exclude_dates,
start = min(.$date),
end = max(.$date),
knots.per.year = 12,
verbose = FALSE)
# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_6x <- deficit_plot(fit_6x,
title = "Deviations from Expected Counts for Pneumonia in MA"
)
deficit_plot_6x
# -- Intervals of inordinate mortality found by the excess model
fit_6x$detected_intervals
## start end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-10-27 2019-12-01 2.532637 2.967751 0.06061043 2046
## 2 2020-04-26 2021-07-04 1.489661 2.618500 0.01756974 12636
## 3 2021-09-05 2022-05-22 1.947460 2.813849 0.02345135 9964
## 4 2022-07-17 2022-08-28 1.694437 2.049678 0.04663405 1597
## 5 2022-11-06 2022-11-13 2.536350 2.922000 0.10416801 683
## 6 2023-01-22 2023-02-19 2.407862 2.978610 0.06651675 1621
## 7 2023-04-23 2023-05-21 2.294970 2.662542 0.06288869 1545
## 8 2024-01-28 2024-02-11 2.619286 2.974941 0.08581986 1058
## expected deficit sd
## 1 2397.5084 -351.5084 48.96436
## 2 22211.3287 -9575.3287 149.03466
## 3 14396.7975 -4432.7975 119.98666
## 4 1931.8138 -334.8138 43.95240
## 5 786.8496 -103.8496 28.05084
## 6 2005.2341 -384.2341 44.77984
## 7 1792.4540 -247.4540 42.33738
## 8 1201.6586 -143.6586 34.66495
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit <- deficit_cumulative(fit_6x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_6 <- cumulative_deficit %>%
ggplot(aes(date)) +
geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
geom_line(aes(y = observed),
color = "white",
size = 1) +
geom_line(aes(y = observed)) +
geom_point(aes(y = observed)) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Date",
y = "Cumulative Deficit Hospitilization-pneumonia",
title = "Cumulative Deficit Hospitilization for pneumonia Cases in MA",
subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_6
library(patchwork)
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")
combined_plot_6_adj <- combined_plot_6 +
theme(legend.position = "bottom",
legend.justification = c(1, 0),
legend.box.just = "right") +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
deficit_plot_6x_adj <- deficit_plot_6x +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_6_adj <- cumulative_deficit_plot_6 +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
final_combined_plot_6 <- combined_plot_6_adj /
deficit_plot_6x_adj /
cumulative_deficit_plot_6_adj +
plot_layout(ncol = 1, heights = c(1, 1, 1))
final_combined_plot_6
ggsave("Pneumonia_combined_figure.png",
plot = final_combined_plot_6,
width = 11,
height = 10.5,
dpi = 300,
bg = "white")
# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1),
lubridate::make_date(2021, 2, 27),
by = "day")
# renal data
renal <- all_cause_state_by_week %>%
filter(category == "renal") %>%
filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
outcome = statecount_numeric)
# excess model with COVID-19 dates excluded
renal_deficit <- compute_expected(renal,
exclude = exclude_dates,
frequency = 52)
## Overall death rate is 2.42.
# Plotting
combined_plot_7 <- ggplot(renal_deficit, aes(x = date)) +
geom_line(aes(y = expected), color = "blue", size = 1.5) + # Expected outcome line
geom_point(aes(y = outcome, color = factor(excluded)), size = 2) + # Observed outcomes
scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
labels = c("Included in Model", "Excluded from Model")) +
labs(title = "Weekly Renal ED Hospitilizations in MA",
x = "Date",
y = "Weekly Hospitilizations",
color = "Data Point Status") +
theme_minimal()
combined_plot_7
# -- Fitting excess model
fit_7x <- renal_deficit %>%
deficit_model(exclude = exclude_dates,
start = min(.$date),
end = max(.$date),
knots.per.year = 12,
verbose = FALSE)
# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_7x <- deficit_plot(fit_7x,
title = "Deviations from Expected Counts for Renal in MA"
)
deficit_plot_7x
# -- Intervals of inordinate mortality found by the excess model
fit_7x$detected_intervals
## start end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-11-10 2019-12-01 1.914332 2.068933 0.06198012 1031
## 2 2020-03-01 2020-07-12 1.798840 2.309561 0.02928593 4844
## 3 2020-09-06 2020-10-04 2.170195 2.326802 0.05879006 1461
## 4 2020-11-01 2021-02-28 1.966115 2.218355 0.03025440 4765
## 5 2021-12-26 2022-01-02 2.213272 2.330685 0.09303280 596
## 6 2022-05-15 2022-05-29 2.364289 2.607537 0.08034593 955
## expected deficit sd
## 1 1114.2633 -83.26328 33.38058
## 2 6219.2915 -1375.29147 78.86248
## 3 1566.4295 -105.42948 39.57814
## 4 5376.3201 -611.32014 73.32339
## 5 627.6176 -31.61763 25.05230
## 6 1053.2544 -98.25443 32.45388
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit <- deficit_cumulative(fit_7x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_7 <- cumulative_deficit %>%
ggplot(aes(date)) +
geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
geom_line(aes(y = observed),
color = "white",
size = 1) +
geom_line(aes(y = observed)) +
geom_point(aes(y = observed)) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Date",
y = "Cumulative Deficit Hospitilization-renal",
title = "Cumulative Deficit Hospitilization for renal Cases in MA",
subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_7
library(patchwork)
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")
combined_plot_7_adj <- combined_plot_7 +
theme(legend.position = "bottom",
legend.justification = c(1, 0),
legend.box.just = "right") +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
deficit_plot_7x_adj <- deficit_plot_7x +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_7_adj <- cumulative_deficit_plot_7 +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
final_combined_plot_7 <- combined_plot_7_adj /
deficit_plot_7x_adj /
cumulative_deficit_plot_7_adj +
plot_layout(ncol = 1, heights = c(1, 1, 1))
final_combined_plot_7
ggsave("Renal_combined_figure.png",
plot = final_combined_plot_7,
width = 11, # good for standard document width
height = 10.5, # leaves space for three stacked plots
dpi = 300,
bg = "white")
# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1),
lubridate::make_date(2021, 6, 12),
by = "day")
# respiratory data
respiratory <- all_cause_state_by_week %>%
filter(category == "respiratory") %>%
filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
outcome = statecount_numeric)
# excess model with COVID-19 dates excluded
respiratory_deficit <- compute_expected(respiratory,
exclude = exclude_dates,
frequency = 52)
## Overall death rate is 3.89.
# Plotting
combined_plot_8 <- ggplot(respiratory_deficit, aes(x = date)) +
geom_line(aes(y = expected), color = "blue", size = 1.5) + # Expected outcome line
geom_point(aes(y = outcome, color = factor(excluded)), size = 2) + # Observed outcomes
scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
labels = c("Included in Model", "Excluded from Model")) +
labs(title = "Weekly Respiratory ED Hospitilizations in MA",
x = "Date",
y = "Weekly Hospitilizations",
color = "Data Point Status") +
theme_minimal()
combined_plot_8
# -- Fitting excess model
fit_8x <- respiratory_deficit %>%
deficit_model(exclude = exclude_dates,
start = min(.$date),
end = max(.$date),
knots.per.year = 12,
verbose = FALSE)
# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_8x <- deficit_plot(fit_8x,
title = "Deviations from Expected Counts for Respiratory in MA"
)
deficit_plot_8x
# -- Intervals of inordinate mortality found by the excess model
fit_8x$detected_intervals
## start end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-07-21 2019-07-21 2.978262 3.366797 0.15813126 401
## 2 2019-10-27 2019-11-24 3.920017 4.324783 0.08015053 2639
## 3 2020-03-08 2021-05-30 2.675694 4.168126 0.02182343 23417
## 4 2021-11-21 2022-03-27 3.632627 4.441859 0.04166917 9293
## 5 2023-04-23 2023-05-07 4.186401 4.696742 0.10783183 1691
## expected deficit sd
## 1 453.3133 -52.31326 21.29115
## 2 2911.4933 -272.49331 53.95826
## 3 36478.3943 -13061.39427 190.99318
## 4 11363.1786 -2070.17861 106.59821
## 5 1897.1403 -206.14030 43.55617
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit <- deficit_cumulative(fit_8x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_8 <- cumulative_deficit %>%
ggplot(aes(date)) +
geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
geom_line(aes(y = observed),
color = "white",
size = 1) +
geom_line(aes(y = observed)) +
geom_point(aes(y = observed)) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Date",
y = "Cumulative Deficit Hospitilization-respiratory",
title = "Cumulative Deficit Hospitilization for respiratory Cases in MA",
subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_8
library(patchwork)
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")
combined_plot_8_adj <- combined_plot_8 +
theme(legend.position = "bottom",
legend.justification = c(1, 0),
legend.box.just = "right") +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
deficit_plot_8x_adj <- deficit_plot_8x +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_8_adj <- cumulative_deficit_plot_8 +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
final_combined_plot_8 <- combined_plot_8_adj /
deficit_plot_8x_adj /
cumulative_deficit_plot_8_adj +
plot_layout(ncol = 1, heights = c(1, 1, 1))
final_combined_plot_8
ggsave("Respiratory_combined_figure.png",
plot = final_combined_plot_8,
width = 11,
height = 10.5,
dpi = 300,
bg = "white")
# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1),
lubridate::make_date(2021, 6, 5),
by = "day")
# other data
other <- all_cause_state_by_week %>%
filter(category == "other") %>%
filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
outcome = statecount_numeric)
# excess model with COVID-19 dates excluded
other_deficit <- compute_expected(other,
exclude = exclude_dates,
frequency = 52)
## Overall death rate is 41.3.
# Plotting
combined_plot_9 <- ggplot(other_deficit, aes(x = date)) +
geom_line(aes(y = expected), color = "blue", size = 1.5) + # Expected outcome line
geom_point(aes(y = outcome, color = factor(excluded)), size = 2) + # Observed outcomes
scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
labels = c("Included in Model", "Excluded from Model")) +
labs(title = "Weekly Other ED Hospitilizations in MA",
x = "Date",
y = "Weekly Hospitilizations",
color = "Data Point Status") +
theme_minimal()
combined_plot_9
# -- Fitting excess model
fit_9x <- other_deficit %>%
deficit_model(exclude = exclude_dates,
start = min(.$date),
end = max(.$date),
knots.per.year = 12,
verbose = FALSE)
# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_9x <- deficit_plot(fit_9x,
title = "Deviations from Expected Counts for Other in MA"
)
deficit_plot_9x
# -- Intervals of inordinate mortality found by the excess model
fit_9x$detected_intervals
## start end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2020-03-01 2021-05-23 35.68635 41.60681 0.06895008 312318
## 2 2021-08-22 2021-10-10 40.63174 42.79761 0.19933061 43766
## 3 2021-11-21 2022-02-13 36.00309 40.26895 0.15167805 63018
## 4 2022-03-06 2022-04-03 41.10001 42.97170 0.25264781 27669
## 5 2022-04-24 2022-05-08 42.59682 44.33594 0.33130392 17206
## expected deficit sd
## 1 364132.31 -51814.3141 603.4338
## 2 46098.94 -2332.9402 214.7066
## 3 70484.75 -7466.7499 265.4896
## 4 28929.04 -1260.0429 170.0854
## 5 17908.48 -702.4761 133.8226
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit <- deficit_cumulative(fit_9x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_9 <- cumulative_deficit %>%
ggplot(aes(date)) +
geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
geom_line(aes(y = observed),
color = "white",
size = 1) +
geom_line(aes(y = observed)) +
geom_point(aes(y = observed)) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Date",
y = "Cumulative Deficit Hospitilization-other",
title = "Cumulative Deficit Hospitilization for other Cases in MA",
subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_9
library(patchwork)
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")
combined_plot_9_adj <- combined_plot_9 +
theme(legend.position = "bottom",
legend.justification = c(1, 0),
legend.box.just = "right") +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
deficit_plot_9x_adj <- deficit_plot_9x +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_9_adj <- cumulative_deficit_plot_9 +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
final_combined_plot_9 <- combined_plot_9_adj /
deficit_plot_9x_adj /
cumulative_deficit_plot_9_adj +
plot_layout(ncol = 1, heights = c(1, 1, 1))
final_combined_plot_9
ggsave("Other_combined_figure.png",
plot = final_combined_plot_9,
width = 11,
height = 10.5,
dpi = 300,
bg = "white")
library(ggplot2)
library(patchwork)
library(cowplot)
library(scales)
# Common scale settings
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(x_limits[1], x_limits[2], by = "1 year")
category_labels <- c("Cardiovascular", "Digestive", "Infectious", "Infectious Respiratory",
"Injury", "Pneumonia", "Renal", "Respiratory", "Other")
# ----- CLEANERS -----
clean_combined <- function(plots, labels) {
mapply(function(p, label) {
p +
labs(title = label, x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y") +
theme_minimal() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5, size = 11),
axis.text.x = element_text(size = 8))
}, plots, labels, SIMPLIFY = FALSE)
}
clean_deficit <- function(plot_list, labels) {
mapply(function(p, label) {
p +
labs(title = label, x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, size = 11),
axis.text.x = element_text(size = 8))
}, plot_list, labels, SIMPLIFY = FALSE)
}
clean_cumulative <- function(plots, labels) {
mapply(function(p, label) {
p +
labs(title = label, x = NULL, y = NULL, subtitle = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y") +
scale_y_continuous(labels = comma) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, size = 11),
axis.text = element_text(size = 8))
}, plots, labels, SIMPLIFY = FALSE)
}
# ----- APPLY CLEANING -----
combined_clean <- clean_combined(list(
combined_plot_1, combined_plot_2, combined_plot_3,
combined_plot_4, combined_plot_5, combined_plot_6,
combined_plot_7, combined_plot_8, combined_plot_9
), category_labels)
deficit_clean <- clean_deficit(list(
deficit_plot_x, deficit_plot_2x, deficit_plot_3x,
deficit_plot_4x, deficit_plot_5x, deficit_plot_6x,
deficit_plot_7x, deficit_plot_8x, deficit_plot_9x
), category_labels)
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_clean <- clean_cumulative(list(
cumulative_deficit_plot_1, cumulative_deficit_plot_2, cumulative_deficit_plot_3,
cumulative_deficit_plot_4, cumulative_deficit_plot_5, cumulative_deficit_plot_6,
cumulative_deficit_plot_7, cumulative_deficit_plot_8, cumulative_deficit_plot_9
), category_labels)
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
# ----- LEGEND FOR COMBINED -----
legend_shared <- get_legend(
combined_plot_1 +
theme(legend.position = "bottom",
legend.title = element_text(size = 10),
legend.text = element_text(size = 9))
)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
## Warning in get_plot_component(plot, "guide-box"): Multiple components found;
## returning the first one. To return all, use `return_all = TRUE`.
# ----- FINAL COMBINED PLOT (balanced, readable) -----
combined_final <- wrap_plots(combined_clean, ncol = 3) +
plot_layout(guides = "collect") &
theme(legend.position = "bottom")
combined_final <- combined_final +
plot_annotation(
title = "Weekly Observed vs Expected Hospitalizations by Cause",
theme = theme(plot.title = element_text(hjust = 0.5, size = 15))
)
ggsave("combined_all_categories.png",
plot = combined_final,
width = 12, height = 10, dpi = 300, bg = "white")
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
deficit_final <- wrap_plots(deficit_clean, ncol = 3) +
plot_annotation(
title = "Deviations from Expected Weekly Admissions",
theme = theme(plot.title = element_text(hjust = 0.5, size = 14))
)
ggsave("deficit_all_categories.png",
plot = deficit_final,
width = 11, height = 9.5, dpi = 300, bg = "white")
cumulative_grid <- wrap_plots(cumulative_clean, ncol = 3)
y_label <- ggdraw() +
draw_label("Cumulative Deficit in Hospitalizations", angle = 90, size = 12)
cumulative_final <- plot_grid(
y_label, cumulative_grid,
rel_widths = c(0.04, 0.96),
nrow = 1
)
cumulative_final <- cumulative_final +
plot_annotation(
title = "Cumulative Deficit in Hospital Admissions by Cause (2020–2024)",
theme = theme(plot.title = element_text(hjust = 0.5, size = 15))
)
ggsave("cumulative_all_categories.png",
plot = cumulative_final,
width = 12, height = 10, dpi = 300, bg = "white")
combined_final
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
deficit_final
cumulative_final
ggsave("deficit_plot_x.png", plot = deficit_plot_x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_2x.png", plot = deficit_plot_2x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_3x.png", plot = deficit_plot_3x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_4x.png", plot = deficit_plot_4x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_5x.png", plot = deficit_plot_5x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_6x.png", plot = deficit_plot_6x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_7x.png", plot = deficit_plot_7x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_8x.png", plot = deficit_plot_8x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_9x.png", plot = deficit_plot_9x, width = 8, height = 6, dpi = 300)
cumulative_deficit <- deficit_cumulative(fit_x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
direct_plotly <- plot_ly(data = cumulative_deficit) %>%
add_ribbons(x = ~date,
ymin = ~observed - 2*sd,
ymax = ~observed + 2*sd,
name = "95% CI",
fillcolor = "rgba(180, 180, 180, 0.5)",
line = list(color = "transparent"),
hoverinfo = "none") %>%
add_lines(x = ~date,
y = ~observed,
name = "Observed",
line = list(color = "black", width = 2),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
add_markers(x = ~date,
y = ~observed,
name = "Data Points",
marker = list(color = "black", size = 4),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
layout(title = list(text = "Cumulative Deficit Hospitalization for Cardio Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
xaxis = list(title = "Date"),
yaxis = list(title = "Cumulative Deficit Hospitalization-Cardio",
tickformat = ","),
hoverlabel = list(bgcolor = "white",
font = list(family = "Arial", size = 12)),
hovermode = "closest")
direct_plotly
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit <- deficit_cumulative(fit_2x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_2 <- cumulative_deficit %>%
ggplot(aes(date)) +
geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
geom_line(aes(y = observed),
color = "white",
size = 1) +
geom_line(aes(y = observed)) +
geom_point(aes(y = observed)) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Date",
y = "Cumulative Deficit Hospitilization-Digestive",
title = "Cumulative Deficit Hospitilization for Digestive Cases in MA",
subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_2
cumulative_deficit <- deficit_cumulative(fit_2x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
direct_plotly <- plot_ly(data = cumulative_deficit) %>%
add_ribbons(x = ~date,
ymin = ~observed - 2*sd,
ymax = ~observed + 2*sd,
name = "95% CI",
fillcolor = "rgba(180, 180, 180, 0.5)",
line = list(color = "transparent"),
hoverinfo = "none") %>%
add_lines(x = ~date,
y = ~observed,
name = "Observed",
line = list(color = "black", width = 2),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
add_markers(x = ~date,
y = ~observed,
name = "Data Points",
marker = list(color = "black", size = 4),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
layout(title = list(text = "Cumulative Deficit Hospitalization for Digestive Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
xaxis = list(title = "Date"),
yaxis = list(title = "Cumulative Deficit Hospitalization-Digestive",
tickformat = ","),
hoverlabel = list(bgcolor = "white",
font = list(family = "Arial", size = 12)),
hovermode = "closest")
direct_plotly
#3
cumulative_deficit <- deficit_cumulative(fit_3x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
direct_plotly <- plot_ly(data = cumulative_deficit) %>%
add_ribbons(x = ~date,
ymin = ~observed - 2*sd,
ymax = ~observed + 2*sd,
name = "95% CI",
fillcolor = "rgba(180, 180, 180, 0.5)",
line = list(color = "transparent"),
hoverinfo = "none") %>%
add_lines(x = ~date,
y = ~observed,
name = "Observed",
line = list(color = "black", width = 2),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
add_markers(x = ~date,
y = ~observed,
name = "Data Points",
marker = list(color = "black", size = 4),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
layout(title = list(text = "Cumulative Deficit Hospitalization for Infectious Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
xaxis = list(title = "Date"),
yaxis = list(title = "Cumulative Deficit Hospitalization-Infectious",
tickformat = ","),
hoverlabel = list(bgcolor = "white",
font = list(family = "Arial", size = 12)),
hovermode = "closest")
direct_plotly
#4
cumulative_deficit <- deficit_cumulative(fit_4x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
direct_plotly <- plot_ly(data = cumulative_deficit) %>%
add_ribbons(x = ~date,
ymin = ~observed - 2*sd,
ymax = ~observed + 2*sd,
name = "95% CI",
fillcolor = "rgba(180, 180, 180, 0.5)",
line = list(color = "transparent"),
hoverinfo = "none") %>%
add_lines(x = ~date,
y = ~observed,
name = "Observed",
line = list(color = "black", width = 2),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
add_markers(x = ~date,
y = ~observed,
name = "Data Points",
marker = list(color = "black", size = 4),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
layout(title = list(text = "Cumulative Deficit Hospitalization for infectious respiratory Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
xaxis = list(title = "Date"),
yaxis = list(title = "Cumulative Deficit Hospitalization-infectious respiratory",
tickformat = ","),
hoverlabel = list(bgcolor = "white",
font = list(family = "Arial", size = 12)),
hovermode = "closest")
direct_plotly
#5
cumulative_deficit <- deficit_cumulative(fit_5x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
direct_plotly <- plot_ly(data = cumulative_deficit) %>%
add_ribbons(x = ~date,
ymin = ~observed - 2*sd,
ymax = ~observed + 2*sd,
name = "95% CI",
fillcolor = "rgba(180, 180, 180, 0.5)",
line = list(color = "transparent"),
hoverinfo = "none") %>%
add_lines(x = ~date,
y = ~observed,
name = "Observed",
line = list(color = "black", width = 2),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
add_markers(x = ~date,
y = ~observed,
name = "Data Points",
marker = list(color = "black", size = 4),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
layout(title = list(text = "Cumulative Deficit Hospitalization for injury Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
xaxis = list(title = "Date"),
yaxis = list(title = "Cumulative Deficit Hospitalization-injury",
tickformat = ","),
hoverlabel = list(bgcolor = "white",
font = list(family = "Arial", size = 12)),
hovermode = "closest")
direct_plotly
#6
cumulative_deficit <- deficit_cumulative(fit_6x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
direct_plotly <- plot_ly(data = cumulative_deficit) %>%
add_ribbons(x = ~date,
ymin = ~observed - 2*sd,
ymax = ~observed + 2*sd,
name = "95% CI",
fillcolor = "rgba(180, 180, 180, 0.5)",
line = list(color = "transparent"),
hoverinfo = "none") %>%
add_lines(x = ~date,
y = ~observed,
name = "Observed",
line = list(color = "black", width = 2),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
add_markers(x = ~date,
y = ~observed,
name = "Data Points",
marker = list(color = "black", size = 4),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
layout(title = list(text = "Cumulative Deficit Hospitalization for pneumonia Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
xaxis = list(title = "Date"),
yaxis = list(title = "Cumulative Deficit Hospitalization-pneumonia",
tickformat = ","),
hoverlabel = list(bgcolor = "white",
font = list(family = "Arial", size = 12)),
hovermode = "closest")
direct_plotly
#7
cumulative_deficit <- deficit_cumulative(fit_7x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
direct_plotly <- plot_ly(data = cumulative_deficit) %>%
add_ribbons(x = ~date,
ymin = ~observed - 2*sd,
ymax = ~observed + 2*sd,
name = "95% CI",
fillcolor = "rgba(180, 180, 180, 0.5)",
line = list(color = "transparent"),
hoverinfo = "none") %>%
add_lines(x = ~date,
y = ~observed,
name = "Observed",
line = list(color = "black", width = 2),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
add_markers(x = ~date,
y = ~observed,
name = "Data Points",
marker = list(color = "black", size = 4),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
layout(title = list(text = "Cumulative Deficit Hospitalization for renal Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
xaxis = list(title = "Date"),
yaxis = list(title = "Cumulative Deficit Hospitalization-renal",
tickformat = ","),
hoverlabel = list(bgcolor = "white",
font = list(family = "Arial", size = 12)),
hovermode = "closest")
direct_plotly
#8
cumulative_deficit <- deficit_cumulative(fit_8x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
direct_plotly <- plot_ly(data = cumulative_deficit) %>%
add_ribbons(x = ~date,
ymin = ~observed - 2*sd,
ymax = ~observed + 2*sd,
name = "95% CI",
fillcolor = "rgba(180, 180, 180, 0.5)",
line = list(color = "transparent"),
hoverinfo = "none") %>%
add_lines(x = ~date,
y = ~observed,
name = "Observed",
line = list(color = "black", width = 2),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
add_markers(x = ~date,
y = ~observed,
name = "Data Points",
marker = list(color = "black", size = 4),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
layout(title = list(text = "Cumulative Deficit Hospitalization for respiratory Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
xaxis = list(title = "Date"),
yaxis = list(title = "Cumulative Deficit Hospitalization-respiratory",
tickformat = ","),
hoverlabel = list(bgcolor = "white",
font = list(family = "Arial", size = 12)),
hovermode = "closest")
direct_plotly
#9
cumulative_deficit <- deficit_cumulative(fit_9x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
direct_plotly <- plot_ly(data = cumulative_deficit) %>%
add_ribbons(x = ~date,
ymin = ~observed - 2*sd,
ymax = ~observed + 2*sd,
name = "95% CI",
fillcolor = "rgba(180, 180, 180, 0.5)",
line = list(color = "transparent"),
hoverinfo = "none") %>%
add_lines(x = ~date,
y = ~observed,
name = "Observed",
line = list(color = "black", width = 2),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
add_markers(x = ~date,
y = ~observed,
name = "Data Points",
marker = list(color = "black", size = 4),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
layout(title = list(text = "Cumulative Deficit Hospitalization for other Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
xaxis = list(title = "Date"),
yaxis = list(title = "Cumulative Deficit Hospitalization-other",
tickformat = ","),
hoverlabel = list(bgcolor = "white",
font = list(family = "Arial", size = 12)),
hovermode = "closest")
direct_plotly
ggsave("cumulative_deficit_plot_1.png", plot = cumulative_deficit_plot_1, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_2.png", plot = cumulative_deficit_plot_2, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_3.png", plot = cumulative_deficit_plot_3, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_4.png", plot = cumulative_deficit_plot_4, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_5.png", plot = cumulative_deficit_plot_5, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_6.png", plot = cumulative_deficit_plot_6, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_7.png", plot = cumulative_deficit_plot_7, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_8.png", plot = cumulative_deficit_plot_8, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_9.png", plot = cumulative_deficit_plot_9, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_1.png", plot = cumulative_deficit_plot_1, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_2.png", plot = cumulative_deficit_plot_2, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_3.png", plot = cumulative_deficit_plot_3, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_4.png", plot = cumulative_deficit_plot_4, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_5.png", plot = cumulative_deficit_plot_5, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_6.png", plot = cumulative_deficit_plot_6, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_7.png", plot = cumulative_deficit_plot_7, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_8.png", plot = cumulative_deficit_plot_8, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_9.png", plot = cumulative_deficit_plot_9, width = 5, height = 3, dpi = 300)
# new 'category' column added
cardiovascular <- fit_x$detected_intervals %>%
mutate(category = "cardiovascular")
digestive <- fit_2x$detected_intervals %>%
mutate(category = "digestive")
infectious <- fit_3x$detected_intervals %>%
mutate(category = "infectious")
infectious_respiratory <- fit_4x$detected_intervals %>%
mutate(category = "infectious respiratory")
injury <- fit_5x$detected_intervals %>%
mutate(category = "injury")
pneumonia <- fit_6x$detected_intervals %>%
mutate(category = "pneumonia")
renal <- fit_7x$detected_intervals %>%
mutate(category = "renal")
respiratory <- fit_8x$detected_intervals %>%
mutate(category = "respiratory")
other <- fit_9x$detected_intervals %>%
mutate(category = "other")
combined_intervals <- bind_rows(
cardiovascular,
digestive,
infectious,
infectious_respiratory,
injury,
pneumonia,
renal,
respiratory,
other
)
combined_intervals <- combined_intervals %>%
select(category, everything())
combined_intervals <- combined_intervals %>%
mutate(across(c(where(is.numeric)), ~round(., 1)))
kable_table <- combined_intervals %>%
kable(
format = "html",
caption = "Detected Intervals by Disease Category",
align = c("l", rep("c", 9)),
col.names = c(
"Category", "Start Date", "End Date", "Observed Count Rate",
"Expected Count Rate", "SD Count Rate", "Observed Counts",
"Expected Counts", "Deficit", "SD"
),
row.names = FALSE
)
kable_table
| Category | Start Date | End Date | Observed Count Rate | Expected Count Rate | SD Count Rate | Observed Counts | Expected Counts | Deficit | SD |
|---|---|---|---|---|---|---|---|---|---|
| cardiovascular | 2019-09-23 | 2019-12-02 | 8.1 | 8.7 | 0.1 | 12062 | 12840.1 | -778.1 | 113.3 |
| cardiovascular | 2020-01-20 | 2020-02-03 | 8.3 | 8.7 | 0.1 | 3343 | 3515.5 | -172.5 | 59.3 |
| cardiovascular | 2020-03-02 | 2021-02-15 | 7.3 | 8.9 | 0.0 | 49893 | 60938.8 | -11045.8 | 246.9 |
| cardiovascular | 2021-12-27 | 2022-01-31 | 8.3 | 9.0 | 0.1 | 6709 | 7290.1 | -581.1 | 85.4 |
| cardiovascular | 2023-05-01 | 2023-05-01 | 9.7 | 9.6 | 0.3 | 1303 | 1289.0 | 14.0 | 35.9 |
| digestive | 2019-09-15 | 2019-10-13 | 4.0 | 4.3 | 0.1 | 2684 | 2869.1 | -185.1 | 53.6 |
| digestive | 2020-03-01 | 2020-09-06 | 3.5 | 4.4 | 0.0 | 13267 | 16748.9 | -3481.9 | 129.4 |
| digestive | 2020-10-11 | 2021-03-14 | 3.7 | 4.2 | 0.0 | 11504 | 12987.3 | -1483.3 | 114.0 |
| digestive | 2021-12-05 | 2022-01-30 | 3.8 | 4.2 | 0.1 | 4556 | 5073.5 | -517.5 | 71.2 |
| digestive | 2022-05-01 | 2022-05-22 | 4.4 | 4.6 | 0.1 | 2353 | 2502.5 | -149.5 | 50.0 |
| digestive | 2022-09-18 | 2022-10-02 | 4.3 | 4.5 | 0.1 | 1722 | 1835.4 | -113.4 | 42.8 |
| infectious | 2020-05-03 | 2020-08-02 | 2.4 | 2.7 | 0.0 | 4495 | 5124.4 | -629.4 | 71.6 |
| infectious | 2020-10-04 | 2021-05-30 | 2.2 | 2.7 | 0.0 | 10517 | 12848.3 | -2331.3 | 113.4 |
| infectious | 2022-01-23 | 2022-02-27 | 2.5 | 2.8 | 0.1 | 2044 | 2250.7 | -206.7 | 47.4 |
| infectious | 2022-12-25 | 2023-01-22 | 2.6 | 2.8 | 0.1 | 1727 | 1898.7 | -171.7 | 43.6 |
| infectious | 2023-06-11 | 2023-07-02 | 2.5 | 2.8 | 0.1 | 1357 | 1490.1 | -133.1 | 38.6 |
| infectious | 2024-10-20 | 2024-12-15 | 2.5 | 2.8 | 0.0 | 3060 | 3354.2 | -294.2 | 57.9 |
| infectious respiratory | 2019-10-13 | 2019-12-01 | 0.6 | 1.3 | 0.0 | 683 | 1355.1 | -672.1 | 36.8 |
| infectious respiratory | 2020-03-29 | 2021-05-30 | 0.1 | 0.9 | 0.0 | 1036 | 7153.0 | -6117.0 | 84.6 |
| infectious respiratory | 2021-10-31 | 2022-03-20 | 0.5 | 1.3 | 0.0 | 1333 | 3752.6 | -2419.6 | 61.3 |
| infectious respiratory | 2023-01-15 | 2023-02-26 | 0.7 | 1.3 | 0.0 | 616 | 1184.6 | -568.6 | 34.4 |
| infectious respiratory | 2023-04-09 | 2023-04-23 | 0.4 | 0.7 | 0.0 | 170 | 283.3 | -113.3 | 16.8 |
| infectious respiratory | 2024-11-03 | 2024-11-10 | 0.6 | 1.0 | 0.1 | 173 | 270.3 | -97.3 | 16.4 |
| injury | 2019-10-20 | 2019-11-24 | 6.8 | 7.3 | 0.1 | 5458 | 5880.9 | -422.9 | 76.7 |
| injury | 2020-01-12 | 2020-09-13 | 6.3 | 7.6 | 0.0 | 30779 | 36951.6 | -6172.6 | 192.2 |
| injury | 2020-10-04 | 2021-02-07 | 6.6 | 7.5 | 0.1 | 16830 | 19168.3 | -2338.3 | 138.4 |
| injury | 2021-03-07 | 2021-04-11 | 7.1 | 7.6 | 0.1 | 5752 | 6123.8 | -371.8 | 78.3 |
| injury | 2021-05-02 | 2021-05-16 | 7.6 | 7.9 | 0.1 | 3070 | 3208.8 | -138.8 | 56.6 |
| injury | 2021-07-25 | 2021-08-08 | 8.0 | 8.4 | 0.1 | 3223 | 3381.7 | -158.7 | 58.2 |
| injury | 2021-12-19 | 2022-01-16 | 7.1 | 7.7 | 0.1 | 4807 | 5171.7 | -364.7 | 71.9 |
| injury | 2022-07-31 | 2022-08-21 | 8.3 | 8.7 | 0.1 | 4453 | 4668.8 | -215.8 | 68.3 |
| pneumonia | 2019-10-27 | 2019-12-01 | 2.5 | 3.0 | 0.1 | 2046 | 2397.5 | -351.5 | 49.0 |
| pneumonia | 2020-04-26 | 2021-07-04 | 1.5 | 2.6 | 0.0 | 12636 | 22211.3 | -9575.3 | 149.0 |
| pneumonia | 2021-09-05 | 2022-05-22 | 1.9 | 2.8 | 0.0 | 9964 | 14396.8 | -4432.8 | 120.0 |
| pneumonia | 2022-07-17 | 2022-08-28 | 1.7 | 2.0 | 0.0 | 1597 | 1931.8 | -334.8 | 44.0 |
| pneumonia | 2022-11-06 | 2022-11-13 | 2.5 | 2.9 | 0.1 | 683 | 786.8 | -103.8 | 28.1 |
| pneumonia | 2023-01-22 | 2023-02-19 | 2.4 | 3.0 | 0.1 | 1621 | 2005.2 | -384.2 | 44.8 |
| pneumonia | 2023-04-23 | 2023-05-21 | 2.3 | 2.7 | 0.1 | 1545 | 1792.5 | -247.5 | 42.3 |
| pneumonia | 2024-01-28 | 2024-02-11 | 2.6 | 3.0 | 0.1 | 1058 | 1201.7 | -143.7 | 34.7 |
| renal | 2019-11-10 | 2019-12-01 | 1.9 | 2.1 | 0.1 | 1031 | 1114.3 | -83.3 | 33.4 |
| renal | 2020-03-01 | 2020-07-12 | 1.8 | 2.3 | 0.0 | 4844 | 6219.3 | -1375.3 | 78.9 |
| renal | 2020-09-06 | 2020-10-04 | 2.2 | 2.3 | 0.1 | 1461 | 1566.4 | -105.4 | 39.6 |
| renal | 2020-11-01 | 2021-02-28 | 2.0 | 2.2 | 0.0 | 4765 | 5376.3 | -611.3 | 73.3 |
| renal | 2021-12-26 | 2022-01-02 | 2.2 | 2.3 | 0.1 | 596 | 627.6 | -31.6 | 25.1 |
| renal | 2022-05-15 | 2022-05-29 | 2.4 | 2.6 | 0.1 | 955 | 1053.3 | -98.3 | 32.5 |
| respiratory | 2019-07-21 | 2019-07-21 | 3.0 | 3.4 | 0.2 | 401 | 453.3 | -52.3 | 21.3 |
| respiratory | 2019-10-27 | 2019-11-24 | 3.9 | 4.3 | 0.1 | 2639 | 2911.5 | -272.5 | 54.0 |
| respiratory | 2020-03-08 | 2021-05-30 | 2.7 | 4.2 | 0.0 | 23417 | 36478.4 | -13061.4 | 191.0 |
| respiratory | 2021-11-21 | 2022-03-27 | 3.6 | 4.4 | 0.0 | 9293 | 11363.2 | -2070.2 | 106.6 |
| respiratory | 2023-04-23 | 2023-05-07 | 4.2 | 4.7 | 0.1 | 1691 | 1897.1 | -206.1 | 43.6 |
| other | 2020-03-01 | 2021-05-23 | 35.7 | 41.6 | 0.1 | 312318 | 364132.3 | -51814.3 | 603.4 |
| other | 2021-08-22 | 2021-10-10 | 40.6 | 42.8 | 0.2 | 43766 | 46098.9 | -2332.9 | 214.7 |
| other | 2021-11-21 | 2022-02-13 | 36.0 | 40.3 | 0.2 | 63018 | 70484.7 | -7466.7 | 265.5 |
| other | 2022-03-06 | 2022-04-03 | 41.1 | 43.0 | 0.3 | 27669 | 28929.0 | -1260.0 | 170.1 |
| other | 2022-04-24 | 2022-05-08 | 42.6 | 44.3 | 0.3 | 17206 | 17908.5 | -702.5 | 133.8 |
library(kableExtra)
html_table <- combined_intervals %>%
kable("html",
escape = FALSE,
caption = "Detected Intervals by Disease Category",
align = c("l", rep("c", ncol(combined_intervals) - 1)),
col.names = c(
"Category", "Start Date", "End Date", "Observed Count Rate",
"Expected Count Rate", "SD Count Rate", "Observed Counts",
"Expected Counts", "Deficit", "SD"
)) %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE,
position = "center")
html_table
| Category | Start Date | End Date | Observed Count Rate | Expected Count Rate | SD Count Rate | Observed Counts | Expected Counts | Deficit | SD | |
|---|---|---|---|---|---|---|---|---|---|---|
| 1…1 | cardiovascular | 2019-09-23 | 2019-12-02 | 8.1 | 8.7 | 0.1 | 12062 | 12840.1 | -778.1 | 113.3 |
| 2…2 | cardiovascular | 2020-01-20 | 2020-02-03 | 8.3 | 8.7 | 0.1 | 3343 | 3515.5 | -172.5 | 59.3 |
| 3…3 | cardiovascular | 2020-03-02 | 2021-02-15 | 7.3 | 8.9 | 0.0 | 49893 | 60938.8 | -11045.8 | 246.9 |
| 4…4 | cardiovascular | 2021-12-27 | 2022-01-31 | 8.3 | 9.0 | 0.1 | 6709 | 7290.1 | -581.1 | 85.4 |
| 5…5 | cardiovascular | 2023-05-01 | 2023-05-01 | 9.7 | 9.6 | 0.3 | 1303 | 1289.0 | 14.0 | 35.9 |
| 1…6 | digestive | 2019-09-15 | 2019-10-13 | 4.0 | 4.3 | 0.1 | 2684 | 2869.1 | -185.1 | 53.6 |
| 2…7 | digestive | 2020-03-01 | 2020-09-06 | 3.5 | 4.4 | 0.0 | 13267 | 16748.9 | -3481.9 | 129.4 |
| 3…8 | digestive | 2020-10-11 | 2021-03-14 | 3.7 | 4.2 | 0.0 | 11504 | 12987.3 | -1483.3 | 114.0 |
| 4…9 | digestive | 2021-12-05 | 2022-01-30 | 3.8 | 4.2 | 0.1 | 4556 | 5073.5 | -517.5 | 71.2 |
| 5…10 | digestive | 2022-05-01 | 2022-05-22 | 4.4 | 4.6 | 0.1 | 2353 | 2502.5 | -149.5 | 50.0 |
| 6…11 | digestive | 2022-09-18 | 2022-10-02 | 4.3 | 4.5 | 0.1 | 1722 | 1835.4 | -113.4 | 42.8 |
| 1…12 | infectious | 2020-05-03 | 2020-08-02 | 2.4 | 2.7 | 0.0 | 4495 | 5124.4 | -629.4 | 71.6 |
| 2…13 | infectious | 2020-10-04 | 2021-05-30 | 2.2 | 2.7 | 0.0 | 10517 | 12848.3 | -2331.3 | 113.4 |
| 3…14 | infectious | 2022-01-23 | 2022-02-27 | 2.5 | 2.8 | 0.1 | 2044 | 2250.7 | -206.7 | 47.4 |
| 4…15 | infectious | 2022-12-25 | 2023-01-22 | 2.6 | 2.8 | 0.1 | 1727 | 1898.7 | -171.7 | 43.6 |
| 5…16 | infectious | 2023-06-11 | 2023-07-02 | 2.5 | 2.8 | 0.1 | 1357 | 1490.1 | -133.1 | 38.6 |
| 6…17 | infectious | 2024-10-20 | 2024-12-15 | 2.5 | 2.8 | 0.0 | 3060 | 3354.2 | -294.2 | 57.9 |
| 1…18 | infectious respiratory | 2019-10-13 | 2019-12-01 | 0.6 | 1.3 | 0.0 | 683 | 1355.1 | -672.1 | 36.8 |
| 2…19 | infectious respiratory | 2020-03-29 | 2021-05-30 | 0.1 | 0.9 | 0.0 | 1036 | 7153.0 | -6117.0 | 84.6 |
| 3…20 | infectious respiratory | 2021-10-31 | 2022-03-20 | 0.5 | 1.3 | 0.0 | 1333 | 3752.6 | -2419.6 | 61.3 |
| 4…21 | infectious respiratory | 2023-01-15 | 2023-02-26 | 0.7 | 1.3 | 0.0 | 616 | 1184.6 | -568.6 | 34.4 |
| 5…22 | infectious respiratory | 2023-04-09 | 2023-04-23 | 0.4 | 0.7 | 0.0 | 170 | 283.3 | -113.3 | 16.8 |
| 6…23 | infectious respiratory | 2024-11-03 | 2024-11-10 | 0.6 | 1.0 | 0.1 | 173 | 270.3 | -97.3 | 16.4 |
| 1…24 | injury | 2019-10-20 | 2019-11-24 | 6.8 | 7.3 | 0.1 | 5458 | 5880.9 | -422.9 | 76.7 |
| 2…25 | injury | 2020-01-12 | 2020-09-13 | 6.3 | 7.6 | 0.0 | 30779 | 36951.6 | -6172.6 | 192.2 |
| 3…26 | injury | 2020-10-04 | 2021-02-07 | 6.6 | 7.5 | 0.1 | 16830 | 19168.3 | -2338.3 | 138.4 |
| 4…27 | injury | 2021-03-07 | 2021-04-11 | 7.1 | 7.6 | 0.1 | 5752 | 6123.8 | -371.8 | 78.3 |
| 5…28 | injury | 2021-05-02 | 2021-05-16 | 7.6 | 7.9 | 0.1 | 3070 | 3208.8 | -138.8 | 56.6 |
| 6…29 | injury | 2021-07-25 | 2021-08-08 | 8.0 | 8.4 | 0.1 | 3223 | 3381.7 | -158.7 | 58.2 |
| 7…30 | injury | 2021-12-19 | 2022-01-16 | 7.1 | 7.7 | 0.1 | 4807 | 5171.7 | -364.7 | 71.9 |
| 8…31 | injury | 2022-07-31 | 2022-08-21 | 8.3 | 8.7 | 0.1 | 4453 | 4668.8 | -215.8 | 68.3 |
| 1…32 | pneumonia | 2019-10-27 | 2019-12-01 | 2.5 | 3.0 | 0.1 | 2046 | 2397.5 | -351.5 | 49.0 |
| 2…33 | pneumonia | 2020-04-26 | 2021-07-04 | 1.5 | 2.6 | 0.0 | 12636 | 22211.3 | -9575.3 | 149.0 |
| 3…34 | pneumonia | 2021-09-05 | 2022-05-22 | 1.9 | 2.8 | 0.0 | 9964 | 14396.8 | -4432.8 | 120.0 |
| 4…35 | pneumonia | 2022-07-17 | 2022-08-28 | 1.7 | 2.0 | 0.0 | 1597 | 1931.8 | -334.8 | 44.0 |
| 5…36 | pneumonia | 2022-11-06 | 2022-11-13 | 2.5 | 2.9 | 0.1 | 683 | 786.8 | -103.8 | 28.1 |
| 6…37 | pneumonia | 2023-01-22 | 2023-02-19 | 2.4 | 3.0 | 0.1 | 1621 | 2005.2 | -384.2 | 44.8 |
| 7…38 | pneumonia | 2023-04-23 | 2023-05-21 | 2.3 | 2.7 | 0.1 | 1545 | 1792.5 | -247.5 | 42.3 |
| 8…39 | pneumonia | 2024-01-28 | 2024-02-11 | 2.6 | 3.0 | 0.1 | 1058 | 1201.7 | -143.7 | 34.7 |
| 1…40 | renal | 2019-11-10 | 2019-12-01 | 1.9 | 2.1 | 0.1 | 1031 | 1114.3 | -83.3 | 33.4 |
| 2…41 | renal | 2020-03-01 | 2020-07-12 | 1.8 | 2.3 | 0.0 | 4844 | 6219.3 | -1375.3 | 78.9 |
| 3…42 | renal | 2020-09-06 | 2020-10-04 | 2.2 | 2.3 | 0.1 | 1461 | 1566.4 | -105.4 | 39.6 |
| 4…43 | renal | 2020-11-01 | 2021-02-28 | 2.0 | 2.2 | 0.0 | 4765 | 5376.3 | -611.3 | 73.3 |
| 5…44 | renal | 2021-12-26 | 2022-01-02 | 2.2 | 2.3 | 0.1 | 596 | 627.6 | -31.6 | 25.1 |
| 6…45 | renal | 2022-05-15 | 2022-05-29 | 2.4 | 2.6 | 0.1 | 955 | 1053.3 | -98.3 | 32.5 |
| 1…46 | respiratory | 2019-07-21 | 2019-07-21 | 3.0 | 3.4 | 0.2 | 401 | 453.3 | -52.3 | 21.3 |
| 2…47 | respiratory | 2019-10-27 | 2019-11-24 | 3.9 | 4.3 | 0.1 | 2639 | 2911.5 | -272.5 | 54.0 |
| 3…48 | respiratory | 2020-03-08 | 2021-05-30 | 2.7 | 4.2 | 0.0 | 23417 | 36478.4 | -13061.4 | 191.0 |
| 4…49 | respiratory | 2021-11-21 | 2022-03-27 | 3.6 | 4.4 | 0.0 | 9293 | 11363.2 | -2070.2 | 106.6 |
| 5…50 | respiratory | 2023-04-23 | 2023-05-07 | 4.2 | 4.7 | 0.1 | 1691 | 1897.1 | -206.1 | 43.6 |
| 1…51 | other | 2020-03-01 | 2021-05-23 | 35.7 | 41.6 | 0.1 | 312318 | 364132.3 | -51814.3 | 603.4 |
| 2…52 | other | 2021-08-22 | 2021-10-10 | 40.6 | 42.8 | 0.2 | 43766 | 46098.9 | -2332.9 | 214.7 |
| 3…53 | other | 2021-11-21 | 2022-02-13 | 36.0 | 40.3 | 0.2 | 63018 | 70484.7 | -7466.7 | 265.5 |
| 4…54 | other | 2022-03-06 | 2022-04-03 | 41.1 | 43.0 | 0.3 | 27669 | 28929.0 | -1260.0 | 170.1 |
| 5…55 | other | 2022-04-24 | 2022-05-08 | 42.6 | 44.3 | 0.3 | 17206 | 17908.5 | -702.5 | 133.8 |
save_kable(html_table, "combined_intervals_table.html")
webshot("combined_intervals_table.html", "combined_intervals_table.png", vwidth = 1200, vheight = 800)
combined_intervals <- combined_intervals %>%
mutate(across(where(is.numeric), ~ round(., 1)))
table_1 <- combined_intervals %>%
slice(1:5) %>%
kbl(
format = "html",
caption = "Detected Intervals by Category (1 of 2)",
align = c("l", rep("c", 9)),
col.names = c(
"Category", "Start Date", "End Date", "Observed Count Rate",
"Expected Count Rate", "SD Count Rate", "Observed Counts",
"Expected Counts", "Deficit", "SD"
)
) %>%
kable_styling(full_width = FALSE)
save_kable(table_1, file = "detected_intervals_1.png", zoom = 3)
## save_kable will have the best result with magick installed.
table_2 <- combined_intervals %>%
slice(6:9) %>%
kbl(
format = "html",
caption = "Detected Intervals by Category (2 of 2)",
align = c("l", rep("c", 9)),
col.names = c(
"Category", "Start Date", "End Date", "Observed Count Rate",
"Expected Count Rate", "SD Count Rate", "Observed Counts",
"Expected Counts", "Deficit", "SD"
)
) %>%
kable_styling(full_width = FALSE)
save_kable(table_2, file = "detected_intervals_2.png", zoom = 3)
## save_kable will have the best result with magick installed.
table_2
| Category | Start Date | End Date | Observed Count Rate | Expected Count Rate | SD Count Rate | Observed Counts | Expected Counts | Deficit | SD |
|---|---|---|---|---|---|---|---|---|---|
| digestive | 2019-09-15 | 2019-10-13 | 4.0 | 4.3 | 0.1 | 2684 | 2869.1 | -185.1 | 53.6 |
| digestive | 2020-03-01 | 2020-09-06 | 3.5 | 4.4 | 0.0 | 13267 | 16748.9 | -3481.9 | 129.4 |
| digestive | 2020-10-11 | 2021-03-14 | 3.7 | 4.2 | 0.0 | 11504 | 12987.3 | -1483.3 | 114.0 |
| digestive | 2021-12-05 | 2022-01-30 | 3.8 | 4.2 | 0.1 | 4556 | 5073.5 | -517.5 | 71.2 |
library(ggplot2)
library(scales)
total_deficits <- data.frame(
category = c("Cardiovascular", "Digestive", "Infectious", "Infectious Respiratory",
"Injury", "Pneumonia", "Renal", "Respiratory", "Other"),
total_deficit = c(
sum(cumulative_deficit_plot_1$data$observed, na.rm = TRUE),
sum(cumulative_deficit_plot_2$data$observed, na.rm = TRUE),
sum(cumulative_deficit_plot_3$data$observed, na.rm = TRUE),
sum(cumulative_deficit_plot_4$data$observed, na.rm = TRUE),
sum(cumulative_deficit_plot_5$data$observed, na.rm = TRUE),
sum(cumulative_deficit_plot_6$data$observed, na.rm = TRUE),
sum(cumulative_deficit_plot_7$data$observed, na.rm = TRUE),
sum(cumulative_deficit_plot_8$data$observed, na.rm = TRUE),
sum(cumulative_deficit_plot_9$data$observed, na.rm = TRUE)
)
)
bar_plot <- ggplot(total_deficits, aes(x = reorder(category, total_deficit), y = total_deficit)) +
geom_bar(stat = "identity", fill = "steelblue") +
scale_y_continuous(labels = comma) +
labs(title = "Total Cumulative Deficits by Category (2020–2024)",
x = "ICD-10 Category",
y = "Total Cumulative Deficit in Hospitalizations") +
coord_flip() +
theme_minimal(base_size = 12)
bar_plot
ggsave("Total_Cumulative_Deficit_by_Category.png",
plot = bar_plot,
width = 8.5,
height = 6,
dpi = 300,
bg = "white")
cardio_deficit <- tail(deficit_cumulative(fit_x, start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
digestive_deficit <- tail(deficit_cumulative(fit_2x, start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
infectious_deficit <- tail(deficit_cumulative(fit_3x, start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
inf_resp_deficit <- tail(deficit_cumulative(fit_4x, start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
injury_deficit <- tail(deficit_cumulative(fit_5x, start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
pneumonia_deficit <- tail(deficit_cumulative(fit_6x, start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
renal_deficit <- tail(deficit_cumulative(fit_7x, start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
respiratory_deficit<- tail(deficit_cumulative(fit_8x, start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
other_deficit <- tail(deficit_cumulative(fit_9x, start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
deficit_summary <- tibble(
Category = c("Cardiovascular", "Digestive", "Infectious", "Infectious Respiratory",
"Injury", "Pneumonia", "Renal", "Respiratory", "Other"),
Cumulative_Deficit = c(cardio_deficit, digestive_deficit, infectious_deficit,
inf_resp_deficit, injury_deficit, pneumonia_deficit,
renal_deficit, respiratory_deficit, other_deficit)
)
bar_plot2 <- ggplot(deficit_summary, aes(x = reorder(Category, -Cumulative_Deficit),
y = Cumulative_Deficit)) +
geom_bar(stat = "identity", fill = "steelblue") +
geom_text(aes(label = scales::comma(round(Cumulative_Deficit))),
vjust = -0.5, size = 3) +
scale_y_continuous(labels = scales::comma) +
labs(
title = "Cumulative Deficit in Hospitalizations by Category",
subtitle = "Massachusetts, March 2020 – December 2024",
x = "ICD-10 Category",
y = "Cumulative Deficit (Observed - Expected)"
) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
bar_plot2
ggsave("Cumulative Deficit in Hospitalizations by Category.png",
plot = bar_plot2,
width = 8.5,
height = 6,
dpi = 300,
bg = "white")
library(tidyverse)
library(readxl)
library(gt)
library(stringr)
# Load and clean data
icd10categories <- read_excel("icd10categories_updated_321.xlsx") %>%
mutate(icd10code = toupper(trimws(icd10code)))
# Function to generate code ranges and extract sort keys
get_code_range_with_sort <- function(codes) {
codes <- sort(unique(codes))
letter <- substr(codes, 1, 1)
num <- as.integer(str_extract(substr(codes, 2, 4), "\\d+"))
df <- tibble(code = codes, letter = letter, num = num)
ranges <- df %>%
group_by(letter) %>%
summarize(
min_num = min(num),
max_num = max(num),
.groups = "drop"
) %>%
mutate(
range = ifelse(min_num == max_num,
paste0(letter, sprintf("%02d", min_num)),
paste0(letter, sprintf("%02d", min_num), "–", letter, sprintf("%02d", max_num))),
sort_key = min_num + (utf8ToInt(letter) * 1000) # sort letter first, then number
)
list(
range = paste(ranges$range, collapse = ", "),
sort_key = min(ranges$sort_key)
)
}
# Apply function and unpack results
icd10_summary <- icd10categories %>%
group_by(category_0321update) %>%
summarize(
`Number of ICD-10 codes` = n_distinct(icd10code),
.groups = "drop",
temp = list(get_code_range_with_sort(icd10code))
) %>%
mutate(
`ICD-10 code groups` = map_chr(temp, "range"),
sort_key = map_dbl(temp, "sort_key")
) %>%
select(-temp) %>%
rename(category = category_0321update) %>%
arrange(sort_key) # Sort by ICD-10 code order
## Warning: There were 3 warnings in `summarize()`.
## The first warning was:
## ℹ In argument: `sort_key = min_num + (utf8ToInt(letter) * 1000)`.
## ℹ In group 0: .
## Caused by warning:
## ! There was 1 warning in `mutate()`.
## ℹ In argument: `sort_key = min_num + (utf8ToInt(letter) * 1000)`.
## Caused by warning in `utf8ToInt()`:
## ! argument should be a character vector of length 1
## all but the first element will be ignored
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 2 remaining warnings.
# Create GT table
gt_table3 <- icd10_summary %>%
select(-sort_key) %>%
gt() %>%
tab_header(
title = md("**Table. ICD-10 Code Groupings and Counts by Category**")
) %>%
cols_width(
`ICD-10 code groups` ~ px(400)
) %>%
tab_options(table.width = pct(100))
gt_table3
| Table. ICD-10 Code Groupings and Counts by Category | ||
| category | Number of ICD-10 codes | ICD-10 code groups |
|---|---|---|
| infectious | 199 | A00–A99, B01–B99 |
| cardiovascular | 100 | I00–I99 |
| infectious respiratory | 13 | J00–J22 |
| pneumonia | 7 | J12–J18 |
| respiratory | 78 | J19–J99 |
| digestive | 25 | K20–K92 |
| injury | 342 | M48–M97, R45, S00–S99, T00–T88, V00–V98, W00–W88, X00–X99, Y00–Y99 |
| renal | 29 | N00–N28 |
| other | 75 | O00–O99, Z31–Z33 |
| covid | 1 | U07 |
# Save the GT table
gtsave(gt_table3, filename = "icd10_category_groupings_continuous.png")
library(tidyverse)
library(readxl)
library(gt)
library(stringr)
# Load and clean data
icd10categories <- read_excel("icd10categories_updated_321.xlsx") %>%
mutate(icd10code = toupper(trimws(icd10code)))
# Modified condensation function with relaxed grouping logic
condense_icd10_codes <- function(codes, max_gap = 2) {
codes <- sort(unique(codes))
letter_part <- substr(codes, 1, 1)
num_part <- as.integer(str_extract(substr(codes, 2, 4), "\\d+"))
df <- tibble(code = codes, letter = letter_part, num = num_part) %>%
arrange(letter, num)
# Relax grouping: allow small gaps (max_gap)
df <- df %>%
arrange(letter, num) %>%
group_by(letter) %>%
mutate(group = cumsum(c(TRUE, diff(num) > max_gap))) %>%
group_by(letter, group) %>%
summarize(
range = if (n() == 1) first(code) else paste0(first(code), "–", last(code)),
sort_val = first(num) + utf8ToInt(first(letter)) * 1000,
.groups = "drop"
)
list(
range = paste(df$range, collapse = ", "),
sort_key = min(df$sort_val)
)
}
# Summarize and apply function
icd10_summary <- icd10categories %>%
group_by(category_0321update) %>%
summarize(
`Number of ICD-10 codes` = n_distinct(icd10code),
temp = list(condense_icd10_codes(icd10code, max_gap = 2)),
.groups = "drop"
) %>%
mutate(
`ICD-10 code groups` = map_chr(temp, "range"),
sort_key = map_dbl(temp, "sort_key")
) %>%
rename(category = category_0321update) %>%
arrange(sort_key) # Alphabetical ICD-10 ordering
# Create GT table
gt_table <- icd10_summary %>%
select(-sort_key, -temp) %>%
gt() %>%
tab_header(
title = md("**Table 1. ICD-10 Code Groupings and Counts by Category**")
) %>%
cols_width(
`ICD-10 code groups` ~ px(500)
) %>%
tab_options(table.width = pct(100))
gt_table
| Table 1. ICD-10 Code Groupings and Counts by Category | ||
| category | Number of ICD-10 codes | ICD-10 code groups |
|---|---|---|
| infectious | 199 | A00–A99, B01–B99 |
| cardiovascular | 100 | I00–I99 |
| infectious respiratory | 13 | J00–J06, J09–J11, J20–J22 |
| pneumonia | 7 | J12–J18 |
| respiratory | 78 | J19, J23–J99 |
| digestive | 25 | K20–K31, K50–K52, K55–K63, K92 |
| injury | 342 | M48, M80, M84, M96–M97, R45, S00–S99, T00–T88, V00–V05, V09–V15, V18–V29, V40–V49, V53–V74, V79–V83, V86–V98, W00–W01, W04–W22, W25–W31, W34, W39–W40, W44–W46, W49–W57, W61, W67–W69, W74, W86–W88, X00–X02, X08, X14, X19, X30, X58, X74, X78–X80, X83, X95, X99, Y00, Y03–Y04, Y07–Y09, Y21, Y24, Y27–Y32, Y63, Y69–Y73, Y79, Y82–Y84, Y90–Y95, Y99 |
| renal | 29 | N00–N28 |
| other | 75 | O00–O04, O07–O16, O20–O36, O40–O48, O60–O77, O80–O82, O85–O94, O98–O99, Z31–Z33 |
| covid | 1 | U07 |
# Save table
gtsave(gt_table, filename = "icd10_category_groupings.png")